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Abstract 

We present a probabilistic cellular automaton with two absorbing 
states, which can be considered a natural extension of the Domany-Kinzel 
model. Despite its simplicity, it shows a very rich phase diagram, with 
two second-order and one first-order transition lines that meet at a bi- 
critical point. We study the phase transitions and the critical behavior of 
the model using mean field approximations, direct numerical simulations 
and field theory. The second-order critical curves and the kink critical dy- 
namics are found to be in the directed percolation and parity conservation 
universality classes, respectively. The first order phase transition is put in 
evidence by examining the hysteresis cycle. We also study the "chaotic" 
phase, in which two replicas evolving with the same noise diverge, using 
mean field and numerical techniques. Finally, we show how the shape 
of the potential of the field-theoretic formulation of the problem can be 
obtained by direct numerical simulations. 
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1 Introduction 



Probabilistic cellular automata (PCA) have been widely used to model a vari- 
ety of systems with local interactions in physics, chemistry, biology and social 
sciences |[ ||, ||, [|. Moreover, PCA are simple and interesting models that 
can be used to investigate fundamental problems in statistical mechanics. Many 
classical equilibrium spin models can be reformulated as PCA, for example the 
kinetic Ising model with parallel heat-bath dynamics is strictly equivalent to a 
PCA with local parallel dynamics || Q . On the other hand, PCA can be mapped 
to spin models j^] by expressing the transition probabilities as exponentials of a 
local energy. PCA can be used to investigate nonequilibrium phenomena, and in 
particular the problem of phase transitions in the presence of absorbing states. 
An absorbing state is represented by a set of configurations from which the 
system cannot escape, equivalent to an infinite energy well in the language of 
statistical mechanics. A global absorbing state can be originated by one or more 
local transition probabilities which take the value zero or one, corresponding to 
some infinite coupling in the local energy j^]. 

The Domany-Kinzel (DK) model is a boolean PCA on a tilted square lattice 
that has been extensively studied ||, [l(J . Let us denote the two possible states 
of each site with the terms "empty" and "occupied". In this model a site at 
time t is connected to two sites at time 4—1, constituting its neighborhood. 
The control parameters of the model are the local transition probabilities that 
give the probability of having an occupied site at a certain position once given 
the state of its neighborhood. The transition probabilities are symmetric for 
all permutations of the neighborhood, and this property is equivalent to saying 
that they depend on the sum of "occupied" sites in the neighborhood, whence 
the term "totalistic" used to denote this class of automata. 

In the DK model the transition probability from an empty neighborhood to 
an occupied state is zero, thus the empty configuration is an absorbing state. 
For small values of the other transition probabilities, any initial configuration 
will evolve to the absorbing state. For larger values, a phase transition to an 
active phase, represented by an ensemble of partially occupied configurations, is 
found. The order parameter of this transition is the asymptotic average fraction 
of occupied sites, which we call the density. The critical properties of this phase 
transition belong to the directed percolation (DP) universality class (except for 
one extreme point) jy], and the DK model is often considered the prototype of 
such a class. 

The evolution of this kind of models is the discrete equivalent of the tra- 
jectory of a stochastic dynamical system. One can determine the sensitivity 
with respect to a perturbation, by studying the trajectories originating by two 
initially different configurations (replicas) evolving with the same realization of 
the stochasticity, e.g. using the same sequence of random numbers. The order 
parameter here is the asymptotic difference between the two replicas, which we 
call the damage. It turns out that, inside the active phase, there is a "chaotic" 
phase in which the trajectories depend on the initial configurations and the 
damage is different from zero, and a non chaotic one in which all trajectories 
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eventually synchronize with the vanishing of the damage. In simple models like 
the DK one, this transition does not depend on the choice of the initial config- 
urations (provided they are different from the absorbing state) and the initial 
damage Q. 

It has been conjectured that all second-order phase transitions from an "ac- 
tive" phase to a non degenerate, quiescent phase (generally represented by an 
absorbing state) belong to the DP universality class if the order parameter is 
a scalar and there are no extra symmetries or conservation laws . This 

has been verified in a wide class of models, even multi-component, and in the 



presence of several asymmetric absorbing states [16 . Also the damage phase 
transition has a similar structure. Once synchronized, the two replicas cannot 
separate, and thus the synchronized state is absorbing. Indeed, numerical sim- 
ulations shows that it is in the DP universality class fl3]| . Moreover, in the DK 
model, the damage phase transition can be mapped onto the density one |Q. 

On the other hand, some models with conserved quantities Jl7], [l^] or sym- 
metric absorbing states belong to a different universality class called parity 
conservation (PC) or directed Ising (l^, [20| . This universality class appears to 
be less robust since it is strictly related to the symmetry of the absorbing states; 
a slight asymmetry is sufficient to bring the model to the usual DP class |L9| |2(j| . 

An interesting question concerns the simplest, one-dimensional PCA model 
with short range interactions exhibiting a first order phase transition. Dickman 
and Tome pl| , p2[ proposed a contact process with spontaneous annihilation, 
auto catalytic creation by trimers and hopping. They found a first order tran- 
sition for high hopping probability, i.e., in the region more similar to mean field 
(weaker spatial correlations). 

Bassler and Browne discussed a model whose phase diagram also presents 
first and second-order phase transitions |30|. In it, monomers of three different 
chemical species can be adsorbed on a one-dimensional surface and neighboring 
monomers belonging to different species annihilate instantaneously. The control 
parameters of the model are the absorption rates of the monomers. The tran- 
sition from a saturate to a reactive phase belongs to the DP universality class, 
while the transition between two saturated phases is discontinuous. The point 
at which three phase transition lines join does belong to the PC universality 
class. 

Scaling and fluctuations near first order phase transitions are also an inter- 
esting subject of study [^3[ |24|, |5|, [2(|, which can profit from the existence of 
simple models. 

In this paper we study a one-dimensional, one-component, totalistic PCA 
with two absorbing states. It can be considered as a natural extension of the 
DK model to a lattice in which the neighborhood of a site at time t contains 
the site itself and its two nearest neighbors at time t — 1. This space-time 
lattice arises naturally in the discretization of one-dimensional reaction-diffusion 
systems. In our model, the transition probabilities from an empty neighborhood 
is zero, and that from a completely occupied neighborhood is one. The model 
has two absorbing states: the completely empty and the completely occupied 
configurations. The order parameter is again the density; it is zero or one in 
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the two quiescent phases, and assumes other values in the active phase. The 
system presents a line of symmetry in the phase diagram, over which the two 
absorbing phases have the same importance. A more detailed illustration of the 
model can be found in Section || 

This model can arise as a particular case of a nonequilibrium wetting of a 
surface. In this framework, only a single layer of particles can be absorbed on 
the surface. If we assume that particles can be absorbed or desorbed only near 
the boundaries of a patch of already absorbed particles (when the neighborhood 
is not homogeneous), then the completely empty and occupied configurations 
are absorbing states. 

This totalistic PCA can also be interpreted as a simple model of opinion 
formation. It assumes that an individual may change his mind according to 
himself and his two nearest neighbors. The role of social pressure is twofold. If 
there is homogeneity of opinions, individuals cannot disagree (absorbing states), 
otherwise they can agree or disagree with the majority with a certain probability. 

The density phase diagram shows two second order phase transition curves 
separating the quiescent phases from the active one, and a first order transition 
line between the two quiescent phases, as discussed in Section |[ These curves 
meet on the line of symmetry in a bicritical point. We use both mean field 
approximations and direct numerical simulations. The former simple approxi- 
mation gives a qualitatively correct phase diagram. The numerical experiments 
are partially based on the fragment method This is a parallel algorithm 

that implements directly the evolution rule for different values of the control 
parameters on the bits of one or more computer words. 

In Section Q we investigate numerically the second order phase transitions 
and find they belong to the DP universality class. Along the line of symmetry 
of the model the two absorbing phases are equivalent. In Appendix [B| we show 
that on this line one can reformulate the problem in terms of the dynamics 
of kinks between patches of empty and occupied sites. Since the kinks are 
created and annihilated in pairs, the dynamics conserves the initial number of 
kinks modulo two. In this way we can present an exact mapping between a 
model with symmetric absorbing phases and one with parity conservation. We 
find that the critical kink dynamics at the bicritical point belongs to the PC 
universality class. 

In Section^ we study the chaotic phase, using dynamic mean field techniques 
(reported in Appendix |a|) and direct numerical simulations. The location of 
this phase is similar to that of the DK model: it joins the second-order critical 
curves at the boundary of the phase diagram. 

Our model exhibits a first-order phase transition along the line of symmetry 
in the upper part of the phase diagram. A first-order transition is usually 
associated to an hysteresis cycle. It is possible to observe such a phenomena by 
adding a small perturbing field to the absorbing states, as discussed in Section |[ 

The DP universality class is equivalent to the Reggeon field theory ]27j] , which 
in d = corresponds to a quadratic potential with a logarithmic divergence at 
the origin. The Langevin description for systems in the PC class yields a similar 
potential, except for irrelevant terms p8). It has been shown pq] that one can 
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reconstruct the potential from the numerical integration of the Langevin equa- 
tion, which, however, requires special techniques in the presence of absorbing 
states [^9| . In Section ^ we show how the potential is reconstructed from actual 
simulations of a phenomenological model, such as our original cellular automata 
or the kink dynamics. In this way we obtain the shape of the potential for a 
system in the parity conservation universality class. 



We describe here a one-dimensional, totalistic, probabilistic cellular automa- 
ton with three inputs. The state of the model at time t is given by x l — 
(xq, . . . , a4,_x) with x\ S {0, 1}; t = 1,2,... and L is the number of sites. All 
operations on spatial indices are assumed to be modulo L (periodic boundary 
conditions). For simplicity of notation, we write x = x\, x_ — x\_ x , x + — x\ +1 
and x' = We shall indicate by a = x_ + x + x + number of occupied cells 

in the neighborhood. The most general three-input totalistic PCA is defined by 
the quantities p s , which are the conditional probabilities that x' = 1 if a = s. 
The microscopic dynamics of the model is completely specified by 



In this expression R s is a stochastic binary variable that takes the value 1 with 
probability p s and with probability 1 — p s , and 6 is the Kronecker delta. 
In practice, we implement R s by extracting a random number r s uniformly 
distributed between and 1, and setting R s equal to 1 if r s < p s and otherwise. 
Eq. (|l|) implies the use of four random numbers r s for each site. The evolution 
of a single trajectory is not affected by the eventual correlations among the 
r s , since only one 5 aiS is different from zero. This is not true when computing 
the simultaneous evolution of two or more replicas using the same noise. If 
not otherwise stated, we use only one random number for all the r s . More 
discussions about the choice of random numbers can be found in Section || and 
in Appendix [b|. 

Withpo = andp 3 = 1 the model presents two quiescent phases: the phase 
corresponding to the configuration x = (0, . . . , 0) and the phase 1 corresponding 
to the configuration x — (1, . . . , 1). In this case there are two control parameters, 
Pi andp2, and the model is symmetric under the changes p\ 1— P2, P2 — * 1~ Pi 
and x — > x © 1 where © is the exclusive disjunction (or the sum modulo 2). 

3 Phase diagram 

In order to have a qualitative idea of the behavior of the model, we first study 
the mean-field approximation. If c and c' denote the density of occupied sites 
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Figure 1: Mean- field phase diagram for the density c of active sites. The white 
(black) region corresponds to the phase (phase 1). The levels of grey indicate 
different values of the asymptotic density c (active phase), the lightest corre- 
sponds to < c < 1/4, and the next ones to 1/4 < c < 1/2, 1/2 < c < 3/4 and 
3/4 < c < 1. The two quiescent phases coexist in the hatched region. 



at times t and t + 1 respectively, 

c' = 3pic(l-c) 2 + 3p 2 c 2 (l-c)+c 3 . (2) 
This map has three fixed points, Cq, C\, and c 2 given by 

c = 0, ci = — — — , and c 2 = 1. 

1 + dpi — 6p2 

The asymptotic density will assume one of the latter values according to the 
values of the control parameters and the initial state as we show in Fig. |l|. In the 
square 1/3 < p% < 1, < p% < 2/3, the only stable fixed point is c\. Inside this 
square, on the segments pi — 2/3 = m(pi — 1/3) with m < 0, c\ = 1/(1 —to). The 
first fixed point Co is stable when p\ < 1/3 and c 2 is stable when p 2 > 2/3. There 
is a continuous second-order transition from the quiescent phase to the active 
phase on the segment p\ = 1/3, < p 2 <2/3 and another continuous transition 
from the active to the quiescent phase 1 on the segment 1/3 < p\ < l,p 2 = 2/3. 
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Figure 2: Phase diagram for the density of active sites c by numerical exper- 
iments. One run was performed with L = 10000 and T — 10000. The graph 
shows 64 x 64 values of p\ and p2- The color code is the same as in Fig. |. The 
inset represents the density profile along the dashed line. Two critical phase 
transitions are evident. 



In the hatched region of Fig. |l] Cq and are both stable. Their basins of 
attraction are, respectively, the semi-open intervals [0, C\) and (c 1; 1]. Starting 
from a uniformly distributed random value of c, as time t goes to infinity, c tends 
to Co with probability c\ , and to C2 with probability 1 — c%. Since, for pi +p 2 = 1 , 
ci = 1/2, the segment p\ + p 2 = 1, with < p\ < 1/3 and 2/3 < p2 < 1, is 
similar to a first-order transition line between the phase and the phase 1. 

In Fig. ^ we show the phase diagram of the model obtained numerically 
starting from a random initial state with half of the cells occupied. The scenario 
is qualitatively the same as predicted by the mean-field analysis. In the vicinity 
of the point (pi,P2) — (0, 1) we observe a discontinuous transition from c = 
to c = 1. The two second-order phase-transition curves from the active phase 
to the quiescent phases are symmetric, and the critical behavior of the order 
parameter, c for the lower curve and 1 — c for the upper one, is the same. 

Due to the symmetry of the model the two second-order phase transition 
curves meet at a bicritical point (p t , 1 — p t ) where the first-order phase transition 
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line p\ + P2 = 1, Pi < Pt ends. Crossing the second-order phase boundaries 
on a line parallel to the diagonal p\ = P2, the density c exhibits two critical 
transitions, as shown in the inset of Fig. || Approaching the bicritical point the 
critical region becomes smaller, and corrections to scaling increase. Finally, at 
the transition point the two transitions coalesce into a single discontinuous one. 

4 Critical dynamics and universality classes 

We performed standard dynamic Monte Carlo simulations starting from a single 
site in the origin out of the nearest absorbing state, and measured the average 
number of active sites N(t), the survival probability P(t) and the average square 
distance from origin R (t) (averaged over surviving runs) defined as 

k=l i 

p w^E e (E^( fc )V ( 3 ) 

k=l \ i J 

In these expressions, k labels the different runs and K is the total number 
of runs. The quantity u>i is Xi if the nearest absorbing state is and 1 — x- t 
otherwise; 9 is the Heaviside step function, that assumes the value 1 if its 
argument is greater than 0, and the value if it is smaller than 0. 
At the critical point one has 

N{t)~t n , P{t)~r S , R 2 {t)~t z . 

At the transition point {p\ = 0.6625(3), p 2 = 0), we get r\ = 0.308(5), 5 = 
0.160(2) and z = 1.265(5), in agreement with the best known values for the 
directed percolation universality class [ j32f . 

Near the bicritical point, on the line p\ + P2 = 1, the two absorbing states 
have symmetrical weight. We define a kink yi as yi = Xi © Xi+%. For the 
computation of the critical properties of the kink dynamics, one has to replace 
uii with yi in Eq. (|^) . The evolution equation is derived in Appendix In the 
kink dynamics there is only one absorbing state (the empty state) , corresponding 
to one of the two absorbing states or 1. For p\ < p t the asymptotic value of 
the density of kinks is zero and it starts to grow for p\ > p t . In models with 
multiple absorbing states, dynamical exponents may vary with initial conditions. 
Quantities computed only on survival runs (R 2 (t)) appear to be universal, while 
others (namely P(t) and N(t)) are not p3[ . 

We performed dynamic Monte Carlo simulations starting either from one 
and two kinks. In both cases pt = 0.460(2), but the exponents were found to be 
different. Due to the conservation of the number of kinks modulo two, starting 
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Figure 3: Mean-field damage-sprea ding phase diagram. The diagram has been 
obtained numerically iterating Eq. (A.l). The lightest level of grey corresponds 
to < h< 1/8, the next ones to 1/8 < h< 1/4, 1/4 < h< 3/8 and 3/8 < h< 
1/2 respectively. 



from a single site one cannot observe the relaxation to the absorbing state, and 
thus 6 = 0. In this case rj = 0.292(5), z = 1.153(5). On the other hand, starting 
with two neighboring kinks, we find 77 = 0.00(2), i5 = 0.285(5), and z = 1.18(2). 
These results are consistent with those found by other authors [[lTj H |l9) . 

5 The chaotic phase 

Let us now turn to the sensitivity of the model to a variation in the initial 
configuration, i.e. to the study of damage spreading or, equivalently, to the 
location of the chaotic phase. 

Given two replicas x and y, we define the difference w as w = x © y. The 
damage h is defined as the fraction of sites in which w = 1, i.e. as the Hamming 
distance between the configurations x and y. 

The precise location of this phase transition depends on the particular im- 
plementation of the stochasticity. Since the sum of occupied cells in the neigh- 
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Figure 4: Phase diagram for the damage spreading from direct numerical sim- 
ulations. The color code is that of Fig. || Traces of the second order phase 
transitions are present; they join to the first-order (co = 0.5) phase boundary. 



borhood of x is in general different from that of y, the evolution equation (yj) 
for the two replicas uses different random numbers r s . The correlations among 
these random numbers affect the location of the chaotic phase boundary |j4| . 

We limit our investigation to the case of maximal correlations by using just 
one random number per site, i.e. all r s are the same for all s at the same site. 
This gives the smallest possible chaotic region. Notice that we have to extract a 
random number for all sites, even if in one or both replicas have a neighborhood 
configuration for which the evolution rule is deterministic (if a — or a = 3). 

One can write the mean field equation for the damage by taking into account 
all the possible local configurations of two lattices. The evolution equation for 
the damage depends on the correlations among sites but in the simplest case we 
can assume that h(t + 1) depends only on h{t) and c(t), the density of occupied 
sites. In Appendix |A| we find the evolution equation for the damage in the 
mean-field approximation. In Fig. || we show the phase diagram of the chaotic 
phase in this approximation. There is a qualitative agreement with the mean 
field phase diagram found for the DK model |3^| . 
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In Fig. [| we show the phase diagram for the damage found numerically by 
considering the evolution starting from uncorrelated configurations with initial 
density equal to 0.5. The damage region is shown in shades of grey. Outside 
this region there appear small damaged domains on the other phase boundaries. 
This is due either to the divergence of the relaxation time (second-order transi- 
tions) or to the fact that a small difference in the initial configuration can drive 
the system to a different absorbing state (first-order transitions). The chaotic 
domain near the point (pi, P2) = (1, 0) is stable regardless of the initial density. 
On the line P2 = the critical points of the density and the damage coincide at 
Pi- 



6 First-order phase transition and hysteresis cy- 
cle 

First-order phase transitions are usually associated to a hysteresis cycle due to 
the coexistence of two phases. In the absence of absorbing states, the coexistence 
of two stable phases for the same values of the parameters is a transient effect 
in finite systems, due to the presence of fluctuations. To find the hysteresis loop 
we modify the model slightly by putting po = 1 — p^ = e with In this 

way the empty and fully occupied configurations are no longer absorbing. This 
brings the model back into the class of equilibrium models for which there is no 
phase transition in one dimension but metastable states can nevertheless persist 
for long times. The mean field equation for the density c becomes 

c '=e(l_c) 3 + 3pic(l-c) 2 + 
3p 2 c 2 (l-c) + (l-e)c 3 . 

We study the asymptotic density as p\ and P2 move on a line with slope 
1 inside the hatched region of Fig. [|. For pi close to zero, Eq. (^) has only 
one fixed point, which is stable and close to e. As p\ increases adiabatically (by 
taking c at t — equal the previous value of the fixed point) the new asymptotic 
density will still assume this value even when two more fixed points appear one 
of which is unstable and the other stable and close to one. Eventually the first 
fixed point disappears, and the asymptotic density jumps to the stable fixed 
point close to one. Going backwards on the same line, the asymptotic density 
will be close to one until that fixed point disappears and it will jump back to a 
small value close to zero. By symmetry, the hysteresis loop is centered around 
the line P1+P2 = 1 which we identify as a first-order phase transition line inside 
the hatched region. 

The hysteresis region is found by two methods, the dynamical mean field, 



which extends the mean field approximation to blocks of I sites 1 36 , and direct 
numerical experiments. As stated before it is necessary to introduce a small 
perturbation e = p$ = 1 — p^. We consider lines parallel to the diagonal p\ = 
P2 in the parameter space and increase the value of p\ and P2 after a given 
relaxation time t r up to P2 — 1; afterwards the scanning is reverted down to 
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Figure 5: Profile of the hysteresis region for several values of the noise e 
and relaxation time T according to the local structure approximation with I = 
6. The curves represent the intersections of the hysteresis cycle with c = 0.5 
(horizontal dashed line in the inset). The curves join smoothly at p\ = 0, 
P2 = 1 (not represented). Starting from the out most curve, these correspond 
to T = 500, e = 0.0001; T = 1000, e = 0.0001, and T = 500, e = 0.001. The 
dashed lines represent the mean-field hysteresis region. The inset represents the 
cycle along a line parallel to the diagonal pi = p 2 - 

Pi = 0. The hysteresis region for various values of parameters are reported 
Fig. In numerical simulations, one can estimate the size of the hysteresis 
region by starting with configurations and 1, and measuring the size d of the 
region in which the two simulations disagree. 

7 Reconstruction of the potential 

An important point in the study of systems exhibiting absorbing states is the 
formulation of a coarse-grained description using a Langevin equation. It is 
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Figure 6: Reconstruction of potential V(c) for p2 = 0. We performed 10 4 
runs over a system of 500 sites and computed the probability distribution P(n) 
averaging over 500 time steps after discarding 1000 time steps. We limited to 
small times in order to be able to see the divergence at the origin (the absorbing 
state) together with the other local minimum (the active state). 

generally accepted that DP universal behavior is represented by 



where c is the density field, a and b are control parameters and a is a Gaussian 
noise with correlations (ot(x, t)a(x' , t')) = 5 x ,x'&t,t' ■ The diffusion coefficient has 
been absorbed into the parameters a and b and the time scale. This equation can 
be obtained by a mean-field approximation of the evolution equation keeping 
only the relevant terms. The state c(x,t) = is clearly stationary, but its 
absorbing character is given by the balance between fluctuations, which are of 
order of the field itself, and the "potential" part ac — be 2 (See also Ref. |2S|). 

The role of the absorbing state can be illustrated by taking a sequence of 
equilibrium models whose energy landscape exhibits the coexistence of an in- 
finitely deep well (the absorbing state) and another broad local minimum (cor- 
responding to the "active", disordered state), separated by an energy barrier. 
There is no true stationary active state for such a system (with a finite energy 
barrier), since there is always a probability of jumping into the absorbing state. 
However, the system can survive in a metastable active state for time intervals 
of the order of the inverse of the heigth of the energy barrier. The parame- 
ters controlling the height of the energy barrier are the size of the lattice and 
the length of the simulation: the equilibrium systems are two-dimensional with 
asymmetric interactions in the time direction Q . In the limiting case of an infi- 
nite system, the height of the energy barrier is finite below the transition point, 
and infinite above, when the only physically relevant state is the disordered one. 



dc(x, t) 

m 



ac(x, t) — bc 2 (x, t) + V 2 c(x, t) + y c(x, t)a(x, t), 



(5) 
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It is possible to introduce a zero-dimensional approximation to the model by 
averaging over the time and the space, assuming that the system has entered 
the metastable state. In this approximation, the size of the original systems 
enters through the renormalized coefficients a, b, 

— — ' - = ac(x, t) — bc 2 (x, t) + \J c(x, t)a(x, t), 
at 

where also the time scale has been renormalized. 
The associated Fokker-Planck equation is 

8P(c,t) 8 , t 2\ / \ 13 2 , . 

where P(c, t) is the probability of observing a density c at time t. One possible 
solution is a <5-peak centered at the origin, corresponding to the absorbing state. 

By considering only those trajectories that do not enter the absorbing state 
during the observation time, one can impose a detailed balance condition, whose 
effective agreement with the actual probability distribution has to be checked a 
posteriori. 

A stationary distribution P{c) = exp(— V(c)) corresponds to an effective 
potential V(c) of the form 

V(c) = log(c) - 2ac + be 2 . 

Note that this distribution is not normalizable. One can impose a cutoff for low 
c, making P(c) normalizable. For finite systems the only stationary solution is 
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the absorbing state. However, by increasing the size of the system, one approx- 
imates the limit in which the energy barrier is infinitely hight, the absorbing 
state unreachable and P(c) is the observable distribution. 

In order to find the form of the effective potential for spatially extended sys- 
tems, Munoz 28) numerically integrated Eq. (pj), using the procedure described 
by Dickman [29[. It is however possible to obtain the shape of the effective 
potential from the actual simulations, simply by plotting V(c) = — log(P(c)) 
versus c, where c is the density of the configuration. 

In Fig. H we show the profile of the reconstructed potential V for some values 
of p around the critical value or the infinite system p* on the line q = 0. We 
used rather small systems and followed the evolution for a limited amount of 
time in order to balance the weight of the <5-peak with respect to P(c) (which 
is only metastable). For larger systems the absorbing state is not visible above 
the transition and dominates below it. 

On the line q = the model belongs to the DP universality class. One can 
observe that the curve becomes broader in the vicinity of the critical point, 
in correspondence of the divergence of critical fluctuations \ ~ \p ~ Pc\~ 7 j 
7' = 0.54 |32|. By repeating the same type of simulations for the kink dynamics 
(random initial condition), we obtain slightly different curves, as shown in Fig. ^. 
We notice that all curves have roughly the same width. Indeed, the exponent 7' 
for systems in the PC universality class is believed to be exactly p3], as given 
by the scaling relation |32| 7' = dvx — 2f3. Clearly, much more informations can 
be obtained from the knowledge of P(c), either by direct numerical simulations 
or dynamical mean field trough finite scale analysis, as shown for instance in 
Ref. (37 . 



8 Discussions and conclusions 

We have studied a probabilistic cellular automata with two absorbing states 
and two control parameters. This is a simple and natural extension of the 
Domany-Kinzcl (DK) model. Despite its simplicity it has a rich phase diagram 
with two symmetric second-order phase curves that join a first-order line at a 
bicritical point. The phase diagram and the critical properties of the model were 
found using several mean field approximations and numerical simulations. The 
second-order phase transitions belong to the directed percolation universality 
class except for the bicritical point, which belongs to the parity conservation 
(PC) universality class. The first-order phase transition line was put in evidence 
by a modification of the model that allows one to find the hysteresis cycles. 
The model also presents a chaotic phase analogous to the one present in the 
DK model. This phase was studied using direct numerical simulations and 
dynamical mean field. 

On the line of symmetry of the model the relevant behavior is given by 
kink dynamics. We found a closed expression for the kink evolution rule and 
studied its critical properties, which belong to the PC universality class. The 
effective potential governing the coarse-grained evolution for the DP and the PC 
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phase was found through direct simulations, confirming that critical fluctuations 
diverge at most logarithmically in the PC class. 

The phase diagram of our model is qualitatively similar to Bassler and 



Browne's (BB) one 30 . In both models two critical lines in the DP univer- 
sality class meet at a bicritical point in the PC universality class, and give 
origin to a first-order transition line. This suggests that the observed behavior 
has a certain degree of universality. 

An interesting feature of the BB model is that the absorbing states at the 
bicritical point arc indeed symmetric, but the model does not show any con- 
served quantities. We have shown that the bicritical dynamics of our model can 
be exactly formulated either in terms of symmetric states or of kinks dynamics, 
providing an exact correspondence between the presence of conserved quantities 
and the symmetry of absorbing states. 

Furthermore, in order to obtain a qualitatively correct mean-field phase di- 
agram of the BB model, one has to include correlations between triplets, while 
the mean-field phase diagram of our model is already correct at first approxi- 
mation. This suggests that we have described a simpler model, which can be 
used as prototype for multi-critical systems. 
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A Damage spreading in the mean field approx- 
imation 

The minimum damage spreading occurs when the two replicas x and y evolve 
using maximally correlated random numbers, i.e. when all r s in Eq. (|l|) are the 
same. Q]. Let w — x y be the damage at a site i and time t. It is also 
possible to consider w as an independent variable and write y = x © w. We 
denote s = x_ + x + x + , s' = j/_ + y + y+ = (x- u>_) + (x © w) + {x + ® w+) 
and s" = u>_ + w + u>+. The evolution equation for h, the density of damaged 
sites w at time i, is obtained by considering all the local configurations x-xx + 
and W-WW+ of one replica and of the damage 

h'= K(c,s,3)n{h,s",3)\p s - Ps ,\, (A.l) 

X-XX + 

where 

w(a,n,m) =a"(l-a) n " m . 
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In this equation all the sums run from zero to one. The value of c is given by 
Eq. (Q) . The term \p s — p s i | is the probability that R s © R s > is one using only 
one random number for the r s . The argument of the sum is the probability that 



x' =/= y' . It is possible to rewrite Eq.(A.l) in a different form 



h , =W m H m S W™ !V( C > s,m)«(ri,s + s' -2t,m)\p s -p s ,\, (A.2) 



where s and s' are the same as above, and £ is the overlap between x and y, 
i.e. I — (x- A y) + (x Ay) + (x+ A y+) (A is the AND operation). Assuming 
that (£) = if b > a, a < or b < 0, the sum of ( |A.2| ) can run over all positive 
integers. This expression is valid for all totalistic rules with a neighborhood of 
size m (here m = 3). 



The stationary state of Eq. (A.l) (or Eq. (|A.2|)) can be found analytically 



using a symbolic manipulation program. The chaotic transition line is 



with 1/3 < pi < 1,0 <p 2 < 2/3. 



B Kink dynamics 

On the segment p\ +P2 = 1, Pi < Pt the order parameter is the number of kinks. 
The dynamics of the kinks = Xi ©x.; + i (that for the ease of notation we write 
y — x © x+) is obtained by taking the exclusive disjunction of x' = and 
x' + = x\^_\ given by Eq. (Q). In order to obtain a closed expression for the y, a 
little of Boolean algebra is needed. 

The totalistic functions <5 CTiS where s — x^ + x + x + can be expressed in 
terms of the symmetric polynomials £W of degree j [ p8[ . These are 

X— © x © 

X — X © X — © XX-|_ , 
£(1)£(2) = j, ^. 

The totalistic functions are given by 

<5., 2 = e (2) ©e (3) , 

In the evolution equation (|l]), one has R\ = 1 if r% < pi, and = 1 if 
r2 < P2- On the line pi + P2 = 1 (i.e. P2 = 1 — Pi), -R2 takes the value 1 if 
1 — r-2 > Pi- Choosing 1 — r2 = r% (this choice does not affect the dynamics 



t(3) = 
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of a single replica) we have R2 = Ri ® 1 and Eq. (|l|) becomes, after some 
manipulations, 

x' = J R(e (1) ©e (2) )©e (2) , 
where R = R\. One can easily check that 

£W = y_ ©y ©x, 

and 

£ (2) = y-y®x. 
Finally, we obtain the evolution equation for the y 
y' = x 1 © x' + 

= R(y^ © y © y_y) © R+{y ®y+® yy+) © y-y ®yy+®y (B.3) 
= i?(y_ V y) © i?+(y V y+) © y_y © yy+ © y, 

In this expression V denotes the disjunction operation (OR). The sum modulo 
2 (XOR) of all yi over the lattice is invariant with time, since all repeated terms 
cancel out (a © a = 0). Note that the kink dynamics uses correlated noise 
between neighboring sites. 
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